library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(knitr)
library(ggplot2)
library(cluster)
library(factoextra)
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
require(caret)
## Loading required package: caret
## Loading required package: lattice
library(corrplot)
## corrplot 0.92 loaded
library(Metrics)
##
## Attaching package: 'Metrics'
## The following objects are masked from 'package:caret':
##
## precision, recall
library(boot)
##
## Attaching package: 'boot'
## The following object is masked from 'package:lattice':
##
## melanoma
library(hopkins)
library(e1071)
library(mclust)
## Package 'mclust' version 6.0.0
## Type 'citation("mclust")' for citing this R package in publications.
library(fpc)
library(NbClust)
library(GGally)
## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2
library(FeatureImpCluster)
## Loading required package: data.table
##
## Attaching package: 'data.table'
## The following objects are masked from 'package:dplyr':
##
## between, first, last
library(pheatmap)
library(ggplot2)
options(dplyr.summarise.inform = FALSE)
data <- read.csv("train.csv")
grouped_data <- data %>% group_by(MSZoning)
agg_tbl <- grouped_data %>% summarise(median(SalePrice))
grouped_dataset <- as.data.frame(agg_tbl)
colnames(grouped_dataset)[colnames(grouped_dataset) == "MSZoning"] ="Zona"
colnames(grouped_dataset)[colnames(grouped_dataset) == "median(SalePrice)"] ="Promedio de Precio (USD)"
grouped_dataset <- grouped_dataset[order(grouped_dataset$`Promedio de Precio (USD)`, decreasing = TRUE), ]
kable(grouped_dataset, caption = "Promedio de ventas por cada zona")
| Zona | Promedio de Precio (USD) | |
|---|---|---|
| 2 | FV | 205950 |
| 4 | RL | 174000 |
| 3 | RH | 136500 |
| 5 | RM | 120500 |
| 1 | C (all) | 74700 |
Un 50% de las casas en la zona
Floating Village Residential se venden por arriba de los
205940.00 USD. Posteriormente sigue
Residential Low Density.
grouped_area_data <- data %>% group_by(Neighborhood)
area_with_median <- grouped_area_data %>% summarise(median(LotArea))
grouped_dataset_lot <- as.data.frame(area_with_median)
colnames(grouped_dataset_lot)[colnames(grouped_dataset_lot) == "Neighborhood"] ="Vecindario"
colnames(grouped_dataset_lot)[colnames(grouped_dataset_lot) == "median(LotArea)"] ="Area del terreno ft^2"
grouped_dataset_lot <- grouped_dataset_lot[order(grouped_dataset_lot$`Area del terreno ft^2`, decreasing = TRUE), ]
kable(grouped_dataset_lot[0:5,], caption = "Mediana de area por cada vecindario")
| Vecindario | Area del terreno ft^2 | |
|---|---|---|
| 5 | ClearCr | 17575.0 |
| 24 | Timber | 12781.5 |
| 16 | NoRidge | 12090.0 |
| 17 | NridgHt | 11694.0 |
| 7 | Crawfor | 11500.0 |
ggplot(grouped_dataset_lot[0:5,], aes(x=Vecindario, y=`Area del terreno ft^2`)) +
geom_bar(stat="identity", fill="steelblue")
La información nos demuestra que Clear Creek es el vecindario que tiene los terrenos con mayor área. Un 50% de los terrenos tiene un área mayor a 17575.0 pies cuadrados.
data_with_remodelation <- data[data$YearBuilt != data$YearRemodAdd,]
En el dataset hay 696 registros de casas que fueron remodeladas. A continuación agrupamos por año la cantidad de remodelaciones.
grouped_remodelation_data <- data_with_remodelation %>% group_by(YearRemodAdd)
year_with_count <- grouped_remodelation_data %>% summarise(total_count=n(),
.groups = 'drop')
grouped_remodelation <- as.data.frame(year_with_count)
colnames(grouped_remodelation)[colnames(grouped_remodelation) == "YearRemoAdd"] ="Año"
colnames(grouped_remodelation)[colnames(grouped_remodelation) == "total_count"] ="Cantidad de casas remodeladas"
grouped_remodelation <- grouped_remodelation[order(grouped_remodelation$`Cantidad de casas remodeladas`, decreasing = TRUE), ]
kable(grouped_remodelation[0:5,], caption = "Numero de casas Remodeladas por año")
| YearRemodAdd | Cantidad de casas remodeladas | |
|---|---|---|
| 1 | 1950 | 166 |
| 44 | 2006 | 51 |
| 45 | 2007 | 42 |
| 43 | 2005 | 40 |
| 38 | 2000 | 33 |
Los 5 años con más casas remodeladas son 1950, 2006, 2007, 2005 y 2000 respectivamente. Es importante notar que en 1950 hay mucha más diferencia que con los demás años.
most_remodelated_zone_1950 <- tail(names(sort(table(data_with_remodelation[data_with_remodelation$YearRemodAdd == "1950",]$MSZoning))), 1)
most_remodelated_zone_2006 <- tail(names(sort(table(data_with_remodelation[data_with_remodelation$YearRemodAdd == "2006",]$MSZoning))), 1)
most_remodelated_zone_2007 <- tail(names(sort(table(data_with_remodelation[data_with_remodelation$YearRemodAdd == "2007",]$MSZoning))), 1)
La mayoría de las remodelaciones pertenecen a la zona
RL - Residencial de baja densidad.
grouped_foundation_data <- data %>% group_by(Foundation)
foundation_with_count <- grouped_foundation_data %>% summarise(total_count=n(),
.groups = 'drop')
grouped_foundation <- as.data.frame(foundation_with_count)
colnames(grouped_foundation)[colnames(grouped_foundation) == "Foundation"] ="Material para la base"
colnames(grouped_foundation)[colnames(grouped_foundation) == "total_count"] ="Cantidad de casas"
grouped_foundation <- grouped_foundation[order(grouped_foundation$`Cantidad de casas`, decreasing = TRUE), ]
kable(grouped_foundation[0:5,], caption = "Numero de casas por cada material de base")
| Material para la base | Cantidad de casas | |
|---|---|---|
| 3 | PConc | 647 |
| 2 | CBlock | 634 |
| 1 | BrkTil | 146 |
| 4 | Slab | 24 |
| 5 | Stone | 6 |
La mayoría de casas están hechas de concreto. Puede ser que un tipo menos común, pero más caro, como el ladrillo signifique un precio más alto.
grouped_type_data <- data %>% group_by(BldgType)
type_count <- grouped_type_data %>% summarise(total_count=n(), .groups = 'drop')
grouped_type <- as.data.frame(type_count)
colnames(grouped_type)[colnames(grouped_type) == "BldgType"] ="Tipo de casa"
colnames(grouped_type)[colnames(grouped_type) == "total_count"] ="Cantidad"
grouped_type <- grouped_type[order(grouped_type$`Cantidad`, decreasing = TRUE), ]
kable(grouped_type, caption = "Conteo de tipo de casas")
| Tipo de casa | Cantidad | |
|---|---|---|
| 1 | 1Fam | 1220 |
| 5 | TwnhsE | 114 |
| 3 | Duplex | 52 |
| 4 | Twnhs | 43 |
| 2 | 2fmCon | 31 |
Los tipo de casas más comunes son de 1 familia y de tipo town house.
grouped_data <- data %>% group_by(Neighborhood) %>% summarise(minPrice = min(SalePrice), maxPrice = max(SalePrice))
grouped <- grouped_data %>% mutate(priceRange = maxPrice - minPrice) %>% arrange(desc(priceRange))
colnames(grouped)[colnames(grouped) == "Neighborhood"] = "Vecindario"
colnames(grouped)[colnames(grouped) == "minPrice"] = "Precio Mínimo (USD)"
colnames(grouped)[colnames(grouped) == "maxPrice"] = "Precio Máximo (USD)"
colnames(grouped)[colnames(grouped) == "priceRange"] = "Rango de Precios (USD)"
kable(head(grouped, 3), caption = "Top 3 - Vecindarios con mayor rango de precios")
| Vecindario | Precio Mínimo (USD) | Precio Máximo (USD) | Rango de Precios (USD) |
|---|---|---|---|
| NoRidge | 190000 | 755000 | 565000 |
| NridgHt | 154000 | 611657 | 457657 |
| OldTown | 37900 | 475000 | 437100 |
Los vecindarios con mayor rango de precios son NoRidge, NridgHt y OldTown. Se puede observar que el rango de precios de NridgHt y OldTown es similar, sin embargo es necesario tomar en cuenta que el precio mínimo de OldTown es mucho menor al de NridgHt.
zone_counts <- data %>% count(MSZoning) %>% mutate(Proporcion = prop.table(n) * 100) %>%
rename(Zona = MSZoning, Cantidad = n) %>% mutate(Proporcion = paste0(round(Proporcion, 2), "%")) %>%
arrange(desc(Cantidad))
kable(zone_counts, caption = "Proporción de casas por zona")
| Zona | Cantidad | Proporcion |
|---|---|---|
| RL | 1151 | 78.84% |
| RM | 218 | 14.93% |
| FV | 65 | 4.45% |
| RH | 16 | 1.1% |
| C (all) | 10 | 0.68% |
Se puede observar que la zona en donde hay mayor proporción de casas es la zona RL, seguido de RM tomando en cuenta que la proporción de casas en esta zona es considerablemente menor así como en las zonas FV, RH y C(all).
price_mean_up <- mean(data$SalePrice)
houses_mean_above <- subset(data, SalePrice > price_mean_up)
ggplot(data = houses_mean_above, aes(x = MSZoning)) +geom_bar()
ggplot(data = houses_mean_above, aes(x = BldgType)) +geom_bar()
ggplot(data = houses_mean_above, aes(x = HouseStyle)) +geom_bar()
ggplot(data = houses_mean_above, aes(x = GarageType)) +geom_bar()
Las características más comunes son las siguientes: 1) Zona donde se ubica las casas: RL (Casa residencial de baja densidad) 2) Tipo de vivienda: 1Fam (Unifamiliar) 3) Estilo de la casa: 1Story (Un nivel) 4) Tipo de Garaje: Attchd (Adjuntada a la casa)
print("Counting Pools And Garages")
## [1] "Counting Pools And Garages"
nrow(data) - table(data$PoolArea > 0)
##
## FALSE TRUE
## 7 1453
nrow(data) - nrow(data[data$GarageType=="NA",])
## [1] 1379
print("")
## [1] ""
print("Relation with the seal price")
## [1] "Relation with the seal price"
tapply(data$SalePrice, data$PoolQC > 0, mean)
## TRUE
## 288138.6
tapply(data$SalePrice, data$GarageType=="NA", mean)
## FALSE
## 185479.5
1453 casas tienen psicina 1379 casas cuentan con garaje
Como se puede ver si influyen la relación con el precio, las piscinas suelen elevar el precio de las casas
ggplot(data = data, aes(x = factor(YrSold))) +
geom_bar(fill = "green") +
ggtitle("Ventas de casas por año") +
xlab("Año") +
ylab("Cantidad de ventas")
ggplot(data = data, aes(x = YrSold, y = SalePrice)) +
geom_line(color = "green") +
ggtitle("Tendencia del precio de venta") +
xlab("Año") +
ylab("Precio de venta promedio")
Como se puede ver en las dos gráficas presentadas desde el año 2006 al 2009 existio una tendencia donde subía y bajan las ventas de las casas. Pero desde 2010 este fue en bajando.
columns_used <- c()
neighborhoodNames <- c("NoRidge", "NridgHt", "StoneBr", "Timber", "Veenker", "Somerst", "ClearCr", "Crawfor", "CollgCr", "Blmngtn", "Gilbert", "NWAmes", "SawyerW", "Mitchel", "NAmes", "NPkVill", "SWISU", "Blueste", "Sawyer", "OldTown", "Edwards", "BrkSide", "BrDale", "IDOTRR", "MeadowV")
for(n in 1:length(neighborhoodNames)) {
# Variable minuscula para nuestro uso.
data$neighborhood[data$Neighborhood == neighborhoodNames[n]] <- n
}
columns_used <- append(columns_used, "neighborhood")
hs <- c("1Story", "2Story", "1.5Fin", "SLvl", "SFoyer")
for(n in 1:length(hs)) {
# Variable minuscula para nuestro uso.
data$houseStyle[data$HouseStyle == hs[n]] <- n
}
columns_used <- append(columns_used, "houseStyle")
data$houseZone[data$MSZoning == "A"] <- 1
data$houseZone[data$MSZoning == "C"] <- 2
data$houseZone[data$MSZoning == "FV"] <- 3
data$houseZone[data$MSZoning == "I"] <- 4
data$houseZone[data$MSZoning == "RH"] <- 5
data$houseZone[data$MSZoning == "RL"] <- 6
data$houseZone[data$MSZoning == "RP"] <- 7
data$houseZone[data$MSZoning == "RM"] <- 8
columns_used <- append(columns_used, "houseZone")
data$houseUtilities[data$Utilities == "AllPub"] <- 1
data$houseUtilities[data$Utilities == "NoSewr"] <- 2
data$houseUtilities[data$Utilities == "NoSeWa"] <- 3
data$houseUtilities[data$Utilities == "ELO"] <- 4
columns_used <- append(columns_used, "houseUtilities")
data$roadAccess[data$Condition1 == "Artery"] <- 1
data$roadAccess[data$Condition1 == "Feedr"] <- 2
data$roadAccess[data$Condition1 == "Norm"] <- 3
data$roadAccess[data$Condition1 == "RRNn"] <- 4
data$roadAccess[data$Condition1 == "RRAn"] <- 5
data$roadAccess[data$Condition1 == "PosN"] <- 6
data$roadAccess[data$Condition1 == "PosA"] <- 7
data$roadAccess[data$Condition1 == "RRNe"] <- 8
data$roadAccess[data$Condition1 == "RRAe"] <- 9
columns_used <- append(columns_used, "roadAccess")
data$remodelated[data$YearBuilt != data$YearRemodAdd] <- 1
data$remodelated[data$YearBuilt == data$YearRemodAdd] <- 0
columns_used <- append(columns_used, "remodelated")
data$roofStyle[data$RoofStyle == "Flat"] <- 1
data$roofStyle[data$RoofStyle == "Gable"] <- 2
data$roofStyle[data$RoofStyle == "Gambrel"] <- 3
data$roofStyle[data$RoofStyle == "Hip"] <- 4
data$roofStyle[data$RoofStyle == "Mansard"] <- 5
data$roofStyle[data$RoofStyle == "Shed"] <- 6
columns_used <- append(columns_used, "roofStyle")
data$roofMaterial[data$RoofMatl == "ClyTile"] <- 1
data$roofMaterial[data$RoofMatl == "CompShg"] <- 2
data$roofMaterial[data$RoofMatl == "Membran"] <- 3
data$roofMaterial[data$RoofMatl == "Metal"] <- 4
data$roofMaterial[data$RoofMatl == "Roll"] <- 5
data$roofMaterial[data$RoofMatl == "Tar&Grv"] <- 6
data$roofMaterial[data$RoofMatl == "WdShake"] <- 7
data$roofMaterial[data$RoofMatl == "WdShngl"] <- 8
columns_used <- append(columns_used, "roofMaterial")
data$overallQuality <- data$OverallQual
columns_used <- append(columns_used, "overallQuality")
data$overallCondition <- data$OverallCond
columns_used <- append(columns_used, "overallCondition")
data$exteriorCondition[data$ExterCond == "Po"] <- 1
data$exteriorCondition[data$ExterCond == "Fa"] <- 2
data$exteriorCondition[data$ExterCond == "TA"] <- 3
data$exteriorCondition[data$ExterCond == "Gd"] <- 4
data$exteriorCondition[data$ExterCond == "Ex"] <- 5
columns_used <- append(columns_used, "exteriorCondition")
data$foundationMaterial[data$Foundation == "BrkTil"] <- 1
data$foundationMaterial[data$Foundation == "CBlock"] <- 2
data$foundationMaterial[data$Foundation == "PConc"] <- 3
data$foundationMaterial[data$Foundation == "Slab"] <- 4
data$foundationMaterial[data$Foundation == "Stone"] <- 5
data$foundationMaterial[data$Foundation == "Wood"] <- 6
columns_used <- append(columns_used, "foundationMaterial")
data$basement[is.na(data$BsmtQual)] <- 0
data$basement[!is.na(data$BsmtQual)] <- 1
columns_used <- append(columns_used, "basement")
data$basementCondition[data$BsmtCond == "Ex"] <- 3
data$basementCondition[data$BsmtCond == "Gd"] <- 2
data$basementCondition[data$BsmtCond != "Ex"] <- 1
data$basementCondition[data$BsmtCond != "Gd"] <- 1
data$basementCondition[is.na(data$BsmtCond)] <- 0
columns_used <- append(columns_used, "basementCondition")
data$fireplace[is.na(data$FireplaceQu)] <- 0
data$fireplace[!is.na(data$FireplaceQu)] <- 1
columns_used <- append(columns_used, "fireplace")
data$garageArea <- data$GarageArea
columns_used <- append(columns_used, "garageArea")
data$pool[is.na(data$PoolQC)] <- 0
data$pool[!is.na(data$PoolQC)] <- 1
columns_used <- append(columns_used, "pool")
data$additionalFeature[is.na(data$MiscFeature)] <- 0
data$additionalFeature[!is.na(data$MiscFeature)] <- 1
columns_used <- append(columns_used, "additionalFeature")
data$livingArea <- data$GrLivArea
columns_used <- append(columns_used, "livingArea")
data$yearBuilt <- data$YearBuilt
columns_used <- append(columns_used, "yearBuilt")
data$salePrice <- data$SalePrice
columns_used <- append(columns_used, "salePrice")
tv <- c("WD", "Oth", "New", "ConLw", "ConLI", "ConLD", "Con", "CWD", "COD")
for(n in 1:length(tv)) {
# Variable minuscula para nuestro uso.
data$saleType[data$SaleType == tv[n]] <- n
}
columns_used <- append(columns_used, "saleType")
msz <- c("FV", "RL", "RH", "RM" , "C (all)")
for(n in 1:length(msz)) {
# Variable minuscula para nuestro uso.
data$mSZoning[data$MSZoning == msz[n]] <- n
}
columns_used <- append(columns_used, "mSZoning")
cleanData <- subset(data, select = columns_used)
Agrupamiento de las columnas más importantes
# Definir las columnas más importantes
important_cols <- c("YrSold", "SalePrice", "YearBuilt", "YearRemodAdd", "PoolArea", "GarageArea", "OverallQual", "OverallCond", "LotFrontage", "LotArea", "GarageCars")
# Normalizar variables numericas
cols_num_norm <- data[,important_cols] <- mutate_if(data[,important_cols], is.numeric, scale)
Resumen de las columnas a utilizar para llevar a cabo el clustering
print(summary(data[,important_cols]))
## YrSold.V1 SalePrice.V1 YearBuilt.V1
## Min. :-1.3671863 Min. :-1.838074 Min. :-3.286697
## 1st Qu.:-0.6142282 1st Qu.:-0.641296 1st Qu.:-0.571727
## Median : 0.1387300 Median :-0.225587 Median : 0.057352
## Mean : 0.0000000 Mean : 0.000000 Mean : 0.000000
## 3rd Qu.: 0.8916881 3rd Qu.: 0.416387 3rd Qu.: 0.951306
## Max. : 1.6446462 Max. : 7.226343 Max. : 1.282400
##
## YearRemodAdd.V1 PoolArea.V1 GarageArea.V1
## Min. :-1.6887898 Min. :-0.068668 Min. :-2.212205
## 1st Qu.:-0.8653621 1st Qu.:-0.068668 1st Qu.:-0.647694
## Median : 0.4424348 Median :-0.068668 Median : 0.032833
## Mean : 0.0000000 Mean : 0.000000 Mean : 0.000000
## 3rd Qu.: 0.9268040 3rd Qu.:-0.068668 3rd Qu.: 0.481841
## Max. : 1.2174256 Max. :18.299910 Max. : 4.420012
##
## OverallQual.V1 OverallCond.V1 LotFrontage.V1 LotArea.V1
## Min. :-3.687150 Min. :-4.111561 Min. :-2.01978 Min. :-0.923413
## 1st Qu.:-0.794879 1st Qu.:-0.517023 1st Qu.:-0.45502 1st Qu.:-0.296889
## Median :-0.071812 Median :-0.517023 Median :-0.04324 Median :-0.104028
## Mean : 0.000000 Mean : 0.000000 Mean : 0.00000 Mean : 0.000000
## 3rd Qu.: 0.651256 3rd Qu.: 0.381612 3rd Qu.: 0.40972 3rd Qu.: 0.108671
## Max. : 2.820459 Max. : 3.077516 Max. :10.00422 Max. :20.511245
## NA's :259
## GarageCars.V1
## Min. :-2.3646297
## 1st Qu.:-1.0265059
## Median : 0.3116179
## Mean : 0.0000000
## 3rd Qu.: 0.3116179
## Max. : 2.9878655
##
set.seed(123)
hopkins(data[, important_cols], m=1400)
## Warning in runif(m, min = colmin[i], max = colmax[i]): NAs produced
## [1] 0.9999868
El resultado es de 0.9999868, el cual indica una tendencia a clustering alta ya que es un valor cerca a 1.
data_dist <- dist(data[, important_cols])
fviz_dist(data_dist, show_labels=F)
Como se puede observar el método de hopkins esta en lo correcto y se puede confiar en el resultado dado de VAT ya que muestra patrones de agrupamiento en general.
cols_num_norm_w <- na.omit(cols_num_norm)
fviz_nbclust(cols_num_norm_w, kmeans, method = "gap")
## Warning: did not converge in 10 iterations
## Warning: did not converge in 10 iterations
El numero optimo de clusters es 4
number_clusters <- 4
km <- kmeans(cols_num_norm_w, centers = number_clusters, iter.max = 100)
km$size
## [1] 167 375 332 327
fviz_cluster(km, cols_num_norm_w)
Como se puede ver los datos si estan agrupados, él único que muestra un mínimo de dispersión es el primer cluster (de color rojo)
hc<-hclust(data_dist, method = "ward.D2")
plot(hc, cex=0.5, axes=FALSE)
rect.hclust(hc,k=number_clusters)
También se puede ver que los clusters estan agrupados. En efecto VAT estaba en lo correcto.
Dividimos el test de datos en un 75% para entrenamiento y un 25% para pruebas. Utilizamos el metodo de bootstraping porque nos asegura que la distribución de la data toma en cuenta la población total. Es decir, no tendremos un sesgo genrado por cómo la información está organizada.
set.seed(5)
expectedResult <- cleanData$salePrice
partition <- createDataPartition(y=expectedResult,
p=.75,
list=F)
trainingSet <- cleanData[partition,]
testingSet <- cleanData[-partition,]
Visualizamos la matriz de correlación para ver tendencias en la data.
correlations <- cor(cleanData[,c("neighborhood", "yearBuilt", "overallCondition", "overallQuality", "houseZone", "houseUtilities")], use="pairwise.complete.obs")
corrplot(correlations, method="circle", type="lower", sig.level = 0.01, insig = "blank")
De esta figura obtenemos las siguientes conclusiones: - El año en el que
se hizo la casa está correlacionado con el vecindario en el que se
encuentra. - La condición de la casa tienen una correlación ligeramente
negativa. Puede ser que mientras el año es menor, la condición actual
reduce. - La zona donde se encuentra la casa también está ligeramente
relacionada a la calidad de los materiales.
Continuamos con un segundo diagrama
correlations <- cor(cleanData[,c("neighborhood", "yearBuilt", "overallCondition", "overallQuality", "pool", "fireplace", "foundationMaterial", "exteriorCondition", "additionalFeature", "livingArea" )], use="pairwise.complete.obs")
corrplot(correlations, method="circle", type="lower", sig.level = 0.01, insig = "blank")
Este diagrama nos da otra información de las variables: - Mientras mejor
sea la condición del exterior, mejor es la condición en general de la
casa. - la pisicna está ligeramente correlacionada con tener una fogata
y con la condición del exterior de la casa. - El material usado para la
base de la casa está relacionado con el año de construcción de esta. -
El area de vivienda de la casa incrementa conforme la calidad de la
construcción lo hace. - Una área de vivienda mayor permite tener más
comodidades, como fogata o una piscina.
Si hacemos una gráfica de puntos tomando en cuenta algunas variables importantes podemos obtener más información:
pairs(~yearBuilt+overallQuality+overallCondition+livingArea,data=cleanData,
main="Matriz de correlación")
- la calidad de las casas ha mejorado con el paso de los años. - El área
de vivienda aumenta ligeramente con el paso del tiempo.
Tomando este análisis en cuenta, se utilizarán las siguientes variables como candidatas a mejores predicciones: - yearBuilt - el año de construcción tiene una alta influencia en la calidad de la casa, vecindario donde se encuentra, condición actual y área de vivienda. - livingArea - el area de vivienda es un factor imporante porque indica cuánto espacio puede ser aprovechado en la casa. No es la misma área que necesita una familia de 2 personas a una de 5. Además, un área mayor permite tener comodidades como piscina o fogata. - overallQuality - La calidad de la casa también es importante para el precio. Esto nos indica que el material usado es resistente. - overallCondition - La condición actual de la casa da una mejor impresión al mostrarla a posibles clientes. - pool - Una piscina aumenta el precio dado que el cliente debe tener presupuesto para mantenerla. - fireplace - en Iowa (lugar del dataset) hay inviernos con nieve, los clientes agradecerán esta comodidad
set.seed(5)
livingArea se toma como la variable independiente y se utiliza para hacer el modelo univariado de regresión lineal para predecir el precio de las casas.
singleVariableModel2 <- lm(salePrice~livingArea, data = trainingSet)
summary(singleVariableModel2)
##
## Call:
## lm(formula = salePrice ~ livingArea, data = trainingSet)
##
## Residuals:
## Min 1Q Median 3Q Max
## -337959 -28898 -713 22328 243922
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 16969.633 5029.813 3.374 0.000767 ***
## livingArea 108.156 3.152 34.311 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 53750 on 1095 degrees of freedom
## Multiple R-squared: 0.5181, Adjusted R-squared: 0.5177
## F-statistic: 1177 on 1 and 1095 DF, p-value: < 2.2e-16
Los coeficientes en este modelo son significativos. Sin embargo \(R^2\) no es considerado alto, por tanto el modelo no logra explicar la variación en salePrice en su mayoria.
Ecuación de regresión: \(salePrice = 108.16livingArea + 1.696963\times 10^{4}\)
ggplot(data = trainingSet, mapping = aes(x = livingArea, y = salePrice)) +
geom_point(color = "green", size = 2) +
geom_smooth(method = "lm", se = TRUE) +
labs(title = "Precio de venta ~ Área de vivienda", x = "Área de vivienda", y = "Precio de venta") +
theme_bw() + theme(plot.title = element_text(hjust = 0.5))
## `geom_smooth()` using formula = 'y ~ x'
Ecuación para predecir el precio de venta para el conjunto de prueba.
predSP2<-predict(singleVariableModel2, newdata = testingSet)
head(predSP2)
## 1 4 5 6 8 15
## 201917.1 202674.1 254697.4 164278.6 243016.5 152489.6
length(predSP2)
## [1] 363
head(singleVariableModel2$residuals)
## 2 3 7 9 10 11
## 28037.00676 13363.06021 106813.44777 -78939.06315 -15454.06173 47.72457
plot(singleVariableModel2)
En el gráfico Residuals vs Fitted se puede ver
que los datos se encuentran alrededor de 0 pero no de forma aleatoria
especificamente ya que la mayor cantidad de puntos se centran en un
lugar en específico.
hist(singleVariableModel2$residuals)
boxplot(singleVariableModel2$residuals)
qqnorm(singleVariableModel2$residuals)
qqline(singleVariableModel2$residuals, col="red")
Según los gráficos se puede observar que en el histograma obtenido no
tiene una forma normal, ya que los datos se encuentran distribuidos
hacia la derecha. Por otra parte en el gráfico q-q los extremos se
alejan de la línea y en la caja de bigotes se observa que hay varios
datos que se encuentran en los extremos.
library(nortest)
lillie.test(singleVariableModel2$residuals)
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: singleVariableModel2$residuals
## D = 0.099646, p-value < 2.2e-16
Haciendo una prueba de normalidad el valor de p es menor que 0.05 por tanto la hipótesis nula de normalidad se rechaza y se afirma que los datos de los residuos no siguen una distribución normal.
predSPP2<-predict(singleVariableModel2, newdata = testingSet[,c(19,21)])
library(caret)
RMSE(predSPP2,testingSet$salePrice)
## [1] 62571.91
plot(testingSet$salePrice,col="blue")
points(predSPP2, col="red")
summary(testingSet$salePrice-predSPP2)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -467188.0 -30742.2 -1108.3 178.9 21154.0 339005.7
En base a los resultados anteriores no se puede afirmar que se tenga una buena distribución ya que el RMSE es bastante alto, los datos en el gráfico se encuentran principalmente en la parte inferior del gráfico y no en diagonal. Además, en el resumen de estadisticas no se encuentran alrededor de 0 y tienen valores muy altos.
Utilizamos todas las variables como una primera aproximación.
allVariablesModel <- lm(salePrice ~ ., data=trainingSet)
summary(allVariablesModel)
##
## Call:
## lm(formula = salePrice ~ ., data = trainingSet)
##
## Residuals:
## Min 1Q Median 3Q Max
## -244657 -18333 -2015 15695 203644
##
## Coefficients: (2 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.440e+05 1.261e+05 -4.315 1.75e-05 ***
## neighborhood -2.628e+03 2.782e+02 -9.449 < 2e-16 ***
## houseStyle -2.144e+03 1.088e+03 -1.971 0.048994 *
## houseZone 7.459e+03 2.336e+03 3.193 0.001449 **
## houseUtilities NA NA NA NA
## roadAccess -1.170e+03 1.185e+03 -0.987 0.323691
## remodelated 7.964e+03 2.297e+03 3.467 0.000548 ***
## roofStyle 6.302e+03 1.252e+03 5.034 5.65e-07 ***
## roofMaterial 2.306e+02 1.878e+03 0.123 0.902286
## overallQuality 1.594e+04 1.275e+03 12.501 < 2e-16 ***
## overallCondition 4.513e+03 1.112e+03 4.060 5.26e-05 ***
## exteriorCondition -1.826e+03 3.056e+03 -0.597 0.550408
## foundationMaterial 1.462e+03 2.112e+03 0.692 0.489044
## basement 1.104e+04 6.783e+03 1.628 0.103876
## basementCondition NA NA NA NA
## fireplace 6.022e+03 2.403e+03 2.506 0.012351 *
## garageArea 5.114e+01 6.185e+00 8.269 4.06e-16 ***
## pool 1.218e+05 1.916e+04 6.355 3.10e-10 ***
## additionalFeature 1.013e+03 5.901e+03 0.172 0.863726
## livingArea 5.190e+01 2.877e+00 18.037 < 2e-16 ***
## yearBuilt 2.454e+02 6.453e+01 3.802 0.000152 ***
## saleType 4.877e+02 6.355e+02 0.767 0.442992
## mSZoning -9.137e+03 3.183e+03 -2.870 0.004182 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 32540 on 1046 degrees of freedom
## (30 observations deleted due to missingness)
## Multiple R-squared: 0.8264, Adjusted R-squared: 0.8231
## F-statistic: 249 on 20 and 1046 DF, p-value: < 2.2e-16
El valor de R cuadrado con este modelo es de 0.8231. Un
valor arriba de 0.7 es considerado bueno para problemas de correlación.
(Fernando, 2021) https://www.investopedia.com/terms/r/r-squared.asp#:~:text=In%20finance%2C%20an%20R%2DSquared,depend%20on%20the%20specific%20analysis.
Prediciendo valores con este modelo
predictionAllVariables <- predict(allVariablesModel, testingSet, type="response")
## Warning in predict.lm(allVariablesModel, testingSet, type = "response"):
## prediction from a rank-deficient fit may be misleading
allVariablesOutput <- cbind(testingSet, predictionAllVariables)
Observando primeros valores de la predicción
head(allVariablesOutput[,c("salePrice", "predictionAllVariables")])
## salePrice predictionAllVariables
## 1 208500 216973.9
## 4 140000 214242.5
## 5 250000 299280.3
## 6 143000 159179.9
## 8 200000 224751.9
## 15 157000 160281.3
plot(allVariablesModel)
La gráfica de Residuales vs Fitted no presenta los valores distribuídos
aleatoriamente de manera horizontal centrados en
y=0. Sin
embargo, tampoco está del todo mal si consideramos el valor de R
cuadrado. La gráfica Q-Q también nos da información útil:
paraece que los valores se aproximan a una línea recta entre los
cuartiles teóricos -2 y 2.
complete_predicted <- na.omit(allVariablesOutput)
rmseAllVariables <- rmse(complete_predicted$salePrice,complete_predicted$prediction)
El valor del RMSE es 49481.26.
Graficamos la matríz de correlación para todas las variables del modelo.
modelCorrelations <- cor(allVariablesOutput, method="pearson")
modelCorrelations[is.na(modelCorrelations)] <- 0
corrplot(modelCorrelations, type="upper", order="hclust", tl.col="black", tl.srt=45)
En el resumen generado arriba se puede observar que las variables más significativas son: neighborhood, remodelated, roofStyle, overallQuality, overallCondition, garageArea, pool, livingArea, yearBuilt.
La gráfica de los valores ajustados vs los residuales nos indica que el modelo puede tener un problema de overfitting. Es decir, no logrará generalizar correctamente los valores del set te prueba.
newMultivariableModel <- lm(salePrice ~ neighborhood + remodelated + roofStyle + overallQuality + overallCondition + garageArea + livingArea + yearBuilt, data=trainingSet)
summary(newMultivariableModel)
##
## Call:
## lm(formula = salePrice ~ neighborhood + remodelated + roofStyle +
## overallQuality + overallCondition + garageArea + livingArea +
## yearBuilt, data = trainingSet)
##
## Residuals:
## Min 1Q Median 3Q Max
## -250231 -18907 -2185 15338 281048
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.372e+05 1.135e+05 -5.612 2.53e-08 ***
## neighborhood -2.728e+03 2.539e+02 -10.743 < 2e-16 ***
## remodelated 7.881e+03 2.311e+03 3.410 0.000673 ***
## roofStyle 7.133e+03 1.243e+03 5.740 1.23e-08 ***
## overallQuality 1.607e+04 1.223e+03 13.145 < 2e-16 ***
## overallCondition 5.112e+03 1.011e+03 5.055 5.06e-07 ***
## garageArea 5.169e+01 6.167e+00 8.382 < 2e-16 ***
## livingArea 5.341e+01 2.617e+00 20.403 < 2e-16 ***
## yearBuilt 3.046e+02 5.663e+01 5.379 9.19e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 33470 on 1088 degrees of freedom
## Multiple R-squared: 0.8143, Adjusted R-squared: 0.813
## F-statistic: 596.4 on 8 and 1088 DF, p-value: < 2.2e-16
plot(newMultivariableModel)
La reducción de variables disminuyó el sobreajuste. Sin embargo, también disminuyó el valor de R cuadrado a 0.81. Hay que tomar en cuenta que esto es para los datos de entrenamiento, así que pondremos a prueba al modelo:
newModelPrediction <- predict(newMultivariableModel, testingSet, type="response")
newModelOutput <- cbind(testingSet, newModelPrediction)
Observando primeros valores de la predicción
head(newModelOutput[,c("salePrice", "newModelPrediction")])
## salePrice newModelPrediction
## 1 208500 220360.8
## 4 140000 209398.8
## 5 250000 298290.8
## 6 143000 157311.9
## 8 200000 225137.8
## 15 157000 154551.2
completeNewModelPredicted <- na.omit(newModelOutput)
rmseNewModel <- rmse(completeNewModelPredicted$salePrice,completeNewModelPredicted$newModelPrediction)
El valor del RMSE es 45207.99. Lo cuál es menor al valor
del modelo anterior (49481.26). Esto es un indicador que el nuevo modelo
es mejor generalizando la información.
svmodel7 <- lm(salePrice~livingArea, data = testingSet)
summary(svmodel7)
##
## Call:
## lm(formula = salePrice ~ livingArea, data = testingSet)
##
## Residuals:
## Min 1Q Median 3Q Max
## -452159 -31302 -2239 20302 341908
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 22812.106 9650.226 2.364 0.0186 *
## livingArea 104.457 5.925 17.629 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 62710 on 361 degrees of freedom
## Multiple R-squared: 0.4626, Adjusted R-squared: 0.4611
## F-statistic: 310.8 on 1 and 361 DF, p-value: < 2.2e-16
Podemos ver que \(R^2\) (0.46) es más bajo que cuando se utilizó el conjunto de entrenamiento (\(R^2\)) por tanto la eficiencia del algoritmo para predecir el precio de las casas es menor y la regresión lineal se acopla de menor forma al conjunto de datos.
allVM8 <- lm(salePrice ~ ., data=testingSet)
summary(allVM8)
##
## Call:
## lm(formula = salePrice ~ ., data = testingSet)
##
## Residuals:
## Min 1Q Median 3Q Max
## -227319 -20567 -3701 16035 256371
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.089e+06 3.052e+05 -3.570 0.000411 ***
## neighborhood -2.951e+03 5.956e+02 -4.955 1.16e-06 ***
## houseStyle -3.980e+03 2.489e+03 -1.599 0.110793
## houseZone 1.043e+04 5.596e+03 1.864 0.063216 .
## houseUtilities -1.325e+04 2.349e+04 -0.564 0.573130
## roadAccess 2.003e+03 2.855e+03 0.702 0.483428
## remodelated 1.531e+04 5.531e+03 2.768 0.005962 **
## roofStyle 7.473e+03 2.927e+03 2.553 0.011134 *
## roofMaterial 1.381e+04 3.430e+03 4.025 7.09e-05 ***
## overallQuality 1.611e+04 3.139e+03 5.130 4.95e-07 ***
## overallCondition 2.801e+03 2.778e+03 1.009 0.313937
## exteriorCondition 5.398e+03 8.289e+03 0.651 0.515399
## foundationMaterial -2.176e+03 4.350e+03 -0.500 0.617261
## basement 7.320e+03 1.896e+04 0.386 0.699684
## basementCondition NA NA NA NA
## fireplace 5.562e+03 5.428e+03 1.025 0.306254
## garageArea 2.862e+01 1.542e+01 1.856 0.064295 .
## pool -8.657e+04 2.407e+04 -3.596 0.000372 ***
## additionalFeature 4.398e+02 1.038e+04 0.042 0.966226
## livingArea 5.401e+01 6.367e+00 8.483 7.53e-16 ***
## yearBuilt 5.036e+02 1.547e+02 3.255 0.001251 **
## saleType -1.773e+03 1.933e+03 -0.917 0.359748
## mSZoning -6.299e+03 7.541e+03 -0.835 0.404122
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 43040 on 329 degrees of freedom
## (12 observations deleted due to missingness)
## Multiple R-squared: 0.7533, Adjusted R-squared: 0.7376
## F-statistic: 47.85 on 21 and 329 DF, p-value: < 2.2e-16
El valor de \(R^2\) con este modelo
es de 0.75 y aunque es menor a cuando se utilzó el conjunto
de entrenamiento (\(R^2\) = 0.8231) Aún
es considerado aceptable para problemas de correlación. (Fernando, 2021)
https://www.investopedia.com/terms/r/r-squared.asp#:~:text=In%20finance%2C%20an%20R%2DSquared,depend%20on%20the%20specific%20analysis.
allVM10 <- lm(salePrice ~ neighborhood + remodelated + roofStyle + overallQuality + overallCondition + garageArea + livingArea + yearBuilt, data=testingSet)
summary(allVM10)
##
## Call:
## lm(formula = salePrice ~ neighborhood + remodelated + roofStyle +
## overallQuality + overallCondition + garageArea + livingArea +
## yearBuilt, data = testingSet)
##
## Residuals:
## Min 1Q Median 3Q Max
## -323417 -21063 -4141 17198 293648
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.713e+05 2.589e+05 -2.980 0.00308 **
## neighborhood -3.080e+03 5.396e+02 -5.708 2.42e-08 ***
## remodelated 1.581e+04 5.598e+03 2.825 0.00500 **
## roofStyle 7.315e+03 2.912e+03 2.512 0.01245 *
## overallQuality 1.908e+04 2.857e+03 6.678 9.38e-11 ***
## overallCondition 5.168e+03 2.434e+03 2.124 0.03438 *
## garageArea 3.123e+01 1.462e+01 2.136 0.03334 *
## livingArea 5.017e+01 5.928e+00 8.463 7.01e-16 ***
## yearBuilt 3.716e+02 1.285e+02 2.892 0.00407 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 45310 on 354 degrees of freedom
## Multiple R-squared: 0.7249, Adjusted R-squared: 0.7186
## F-statistic: 116.6 on 8 and 354 DF, p-value: < 2.2e-16
En este caso a pesar de \(R^2\) (0.72) aún es considerado aceptable para problemas de correlación (Fernando, 2021) su valor fue reducido considerablemente vs si se hubiera hecho con el conjunto de datos de entrenamiento (0.81). https://www.investopedia.com/terms/r/r-squared.asp#:~:text=In%20finance%2C%20an%20R%2DSquared,depend%20on%20the%20specific%20analysis.
Es posible concluir que para los modelos tanto univariables como multivariables la eficiencia del modelo se ve reducida y un menor \(R^2\) nos indica que el modelo que el modelo se ajusta menos a los datos y por tanto hay una menor precisión en la predicción de los datos.
| modelo | rCuadrado |
|---|---|
| Lineal | 0.5200 |
| Multivariable | 0.8231 |
| Multivariable Significativas | 0.8130 |
El valor de R cuadrado pareciera ser mejor en el modelo Multivariable que usa todas las variables numérics. Sin embargo, hay que tomar en cuenta que esta medida es para el set de entrenamiento. Si se observa el análisis de residuales para ese modelo se ve que este está sobreajustado.
En azul tenemos el precio correcto y en rojo los valores predecidos para
ese mismo conjunto de datos. Podemos observar que el modelo no es muy
bueno generalizando precios muy altos. Lo cuál confirma la forma que se
ve en la gráfica de residuales.
Este modelo a pesar de ser parecido al anterior puede generaizar mejor
los valores de prueba. Esto se ve claramente en las observaciones con
precios bajos.
Por otro lado, al comparar el RMSE de los modelos (62571.91, 49481.26 y 45207.99), obtenemos que el tercero tiene un menor valor, es decir, es un indicador que es mejor generalizando los resultados.
En conclusión, el modelo que mejor generaliza los precios de las casas es el modelo multivariable que utiliza únicamente las variables importantes: neighborhood, remodelated, roofStyle, overallQuality, overallCondition, garageArea, pool, livingArea, yearBuilt. Hay que tomar en cuenta que podríamos tener mejores resultados utilizando un modelo polinómico en vez de lineal porque algunas características, como livingArea no crecen de forma lineal.